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Abstract — In this paper, a new method is introduced to bUndly 
estimate the transmit power of multiple signal sources in multi- 
antenna fading channels, when the number of sensing devices and 
the number of available samples are sufficiently large compared 
to the number of sources. Recent advances in the field of large 
dimensional random matrix theory are used that result in a 
simple and computationally efficient consistent estimator of the 
power of each source. A criterion to determine the minimum 
number of sensors and the minimum number of samples required 
to achieve source separation is then introduced. Simulations are 
performed that corroborate the theoretical claims and show that 
the proposed power estimator largely outperforms alternative 
power inference techniques. 

Index Terms — Cognitive radio, G-estimation, power estimation, 
random matrix theory, statistical inference. 



I. Introduction 

At a time when radio resources become scarce, the alterna- 
tive offered by cognitive radios [1] is gaining more and more 
interest. A cognitive {or: flexible) wireless network is a set of 
opportunistic entities, referred to as the secondary network, 
that benefit from unused spectrum resources to establish 
communication while generating little or even no interference 
to the licensed networks, collectively referred to as the primary 
network. This is achieved by letting the secondary devices 
sense the communication channel for the presence of active 
transmissions and exchange the collected information among 
the secondary network, in order to perform optimal decisions 
on the opportunistic communication strategy to apply. The 
difficulty for the secondary network does not lie in the de- 
tection of downlink transmissions from fixed access points to 
licensed mobile users in the primary network, but rather in the 
reliable detection of the uplink transmissions from the mobile 
licensed users to the primary access points. If, in addition to 
detecting active transmissions, the secondary devices can, at 
all time, detect the exact number of primary mobile sources 
and evaluate the power used by every individual source, the 

R. Couillet and M. Debbah ai'e with the Alcatel-Lucent Chair on Flexible 
Radio, SUPELEC, Gif sur Yvette, 91192, Plateau de Moulon, 3, Rue Joliot- 
Curie, France e-mail: {romain. couillet, merouane. debbah}@supelec.fr. M. 
Debbah's work is supported by the European Commission, FP7 Network of 
Excellence in Wireless Communications NEWCOM++ and the French ANR 
Project SESAME. 

J. W. Silverstein is with the Department of Mathematics, North Carolina 
State University, Raleigh, North Carolina 27695-8205, jack@math.ncsu.edu. 
J. W. Silverstein's work is supported by the U.S. Army Research Office, Grant 
W911NF-09-1-0266. 

Z. Bai is with the KLAS MOE & School of Mathematics and 
Statistics, Northeast Normal University, Changchun, Jilin 130024. China, 
baizd@nenu.edu.cn. Z. Bai's work is supported by the NSF China, Grant 
10871036, the NUS, Grant R-155-000-079-1 12 and R-155-000-096-646. 



transmission policy in the secondary network can be accurately 
and dynamically adapted. An example of use is found in the 
recent development of femtocells, i.e., small area cells that op- 
erate indoors by overlaying the spectrum licensed to outdoors 
macrocells. Closed access femtocells have the capabilities to 
self-organize and to dynamically access spectrum resources 
[2] -[3]; specifically, the first requirement of a femtocell is to 
minimally interfere the overlaid licensed macrocell network, 
while simultaneously trying to optimize transmission data 
rates within the femtocell. This requires that the femtocells 
be constantly aware of the outdoor activity of the macrocell 
mobile users. As such, macrocell-femtocell networks are cog- 
nitive wireless networks in which the established macrocell 
network is seen as the primary network, while the femtocell 
network plays the role of the opportunistic secondary network. 
In [4], the achievable rates of a two-tier macrocell-femtocell 
network are derived in the very general case where all entities 
in the networks are embedded with multiple antennas. The 
optimal coverage of the secondary networks is computed under 
several assumptions on the side information available at the 
femtocells. Among these assumptions, [4] supposes that the 
femtocells have perfect knowledge of the distances to the 
macrocell user equipments. This last assumption suggests that 
the femtocells have a means, either global positioning system 
or some sort of detection mechanism, to perfectly evaluate the 
distances to the active primary users. In the present work, we 
address the problem of the estimation of the distance of the 
secondary network to each primary user or, more exactly, the 
problem of the estimation of the individual source transmit 
powers. We provide a framework for the secondary network 
(i) to identify the number of primary sources, (ii) to determine 
the number of transmit antennas for every source and (iii) to 
estimate the transmit power from each source. 

The difficulty of estimating transmit powers lies in the little 
information known a priori by the secondary network: the 
transmitted data and the transmission channels are usually 
inaccessible. This has motivated much work in the direction of 
blind signal source detection methods, based on the Neyman- 
Pearson test in Gaussian channels [5], Rayleigh fading chan- 
nels [6], multiple antenna channels [7] and large dimensional 
multi-antenna channels [8], but these successive works are 
designed to answer a binary hypothesis test on the presence or 
absence of a signal source. Alternatively, in [9], a method is 
derived to separate signal sources and estimate the number 
of those sources. To solve the harder problem of power 
inference, it is necessary to assume that the amount of sensors 
in the secondary network is larger than the number of active 
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sources, e.g., individual secondary users are equipped with 
many antennas, or a large number of secondary users, each 
of them equipped with few antennas, collect their received 
data via a central backbone; this assumption is valid in the 
context of femtocells that can communicate through wired 
private or public networks. The condition on the number of 
sensors allows one to model the multi-dimensional channel H 
from the joint primary sources to the secondary users, the joint 
source transmit data X and the additive received noise W as 
large dimensional random matrices with independent entries 
(no specific matrix size definition is required at this point). 
Denoting P a diagonal matrix whose entries are the source 
powers with multiplicities the number of transmit antennas 
of each user, the power detection problem boils down to 
estimating the entries of P from the sole knowledge of the 
received data matrix Y = HP^X + W, as all system 
dimensions (number of antennas per transmit source, number 
of sensors, number of available samples) are large. If the 
available samples largely outnumber the sensors (of several 
orders of magnitude), and the number of sensors are much 
larger than the number of transmit antennas, the strong law 
of large numbers ensures that the diagonal entries of P can 
be retrieved directly from the eigenvalues of YY'^, and the 
problem is immediately solved. When all dimensions are large 
but are of the same order of magnitude, the law of large 
numbers no longer applies and one has to consider results from 
the theory of large dimensional random matrices, e.g., [10], 
used in the present article to derive the asymptotic eigenvalue 
distribution of YY*^ as a function of P. To this day and to the 
best of our knowledge, no computationally-efficient consistent 
estimator for the entries of P has been proposed.^ Among the 
existing techniques are discretization and convex optimization 
strategies [11], [12], which tend to directly invert the result 
from [10] (although an explicit inverse was not available at 
that time), and moment-based approaches [13], [14], which use 
the empirical moments of the eigenvalue distribution YY*^ to 
infer the entries of P. Some of these moment-based methods 
are computationally cheap, but provide in general consistent 
estimators of the moments of the eigenvalue distribution of P, 
instead of estimators of the sought powers. These techniques 
are therefore expected to perform worse than methods that 
would fully exploit the eigenvalue distribution of YY*^, and 
not only a few moments of the distribution. This problem is 
successfully addressed in [15] for the simpler sample covari- 
ance matrix model Y' = P2X, where strongly consistent 
estimators for the individual entries of P are provided, which 
are based on the full eigenvalue distribution of X'^PX. 

The present work generalizes this result to infer the entries 
of P from the observed matrix Y = HP 2 X + W. The novel 
estimator proposed here is strongly consistent with respect 
to growing number of sensors, sources and samples, has a 
very compact form, is computationally efficient and is shown 
in simulations to largely outperform alternative approaches, 
such as moment-based methods. The estimator is moreover 
robust to small system dimensions. We specifically show 

'an estimator Pi of the i"" entry Pi of P is said to be consistent if Pi — 
Pi almost surely when the relevant system dimensions grow large. 
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Fig. 1. A cognitive radio network 

that, if the number of sensing entities is larger than the 
number of active transmitters in the primary network, it is 
possible to evaluate both the exact number of transmitters 
and their respective transmit powers (and, for that matter, the 
number of transmit antennas per source can also be estimated). 
Otherwise, ambiguous scenarios might arise where multiple 
transmitters may be confused as a single transmitter with 
estimated transmit power the average of the true transmit 
powers of these transmitters. Additionally, we provide an 
expression of the minimum number of sensors required to 
separate transmit sources of similar power. 

The remainder of this paper is structured as follows: Section 
II introduces the system model. In Section III, we study the 
asymptotic spectrum of the eigenvalues of YY'^. In Section 
IV, the novel power estimator is derived. Section V provides 
simulation results. Section VI concludes this work. 

Notations: In the following, boldface lower case symbols 
represent vectors, capital boldface characters denote matrices 
(In is the size-A^ identity matrix). The transpose and Hermi- 
tian transpose operators are denoted (•)^ and {■)^, respectively. 
We denote by C+ the set {z e C,3[z] > 0} and by the 
set {z £ C,Sy[z] < 0}. The left-limit in x of a function / is 
denoted /(x— ). 

II. System Model 

Consider a wireless (primary) network in which K entities 
are transmitting data simultaneously on the same frequency 
resource. Transmitter fc S {!,..., K} has transmission power 
Pk and is equipped with Uk antennas. We denote n = 
X^fcLi "fc '■^^ '■^'■^l number of transmit antennas of the primary 
network. Consider also a secondary network composed of a 
total of N, N > n, sensing devices (they may be N single 
antenna devices or multiple devices embedded with multiple 
antennas whose sum equals A^); we shall refer to the A^ sensors 
collectively as the receiver. This scenario is depicted in Figure 
1. To ensure that every sensor in the secondary network, 
e.g., in a femtocell, roughly captures the same amount of 
energy from a given transmitter, we need to assume that 
the respective transmitter-sensor distances are alike. This is a 
realistic assumption for a in-house femtocell network. Denote 
Hfc G ([^Nxuk jjjg multiple antenna channel matrix between 
transmitter k and the receiver. We assume that the entries 
of \/iVHfc are independent and identically distributed with 
zero mean, unit variance and finite fourth order moment. At 
time instant m, transmitter k emits the multi-antenna signal 
^(m) g with entries assumed to be independent and 
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identically distributed of zero mean, unit variance and finite 
fourth order moment. Assume further that at time instant m the 
receive signal is impaired by additive white noise with entries 
of zero mean, variance ct^ and finite fourth order moment on 
every sensor; we denote crw^"*) e the receive noise vector 
where the entries of w^"*^ have unit variance. At time m, the 
receiver therefore senses the signal y^"*' e defined as 



K 



(m) 



(1) 



fe=l 



Assuming the channel fading coefficients are constant over at 
least M consecutive sampUng periods, by concatenating M 
successive signal realizations into Y 
CNxM^ we have 



7(1) 



K 

E 

fe=i 



(2) 



.(1) 



„(M)i 



where = 

w = [w(i),..:,wW] 

rewritten as 

Y = HP5X + (7W 



for every fc, and 
This can be further 



(3) 



where P e M"^" is diagonal with first ni entries Pi, 
subsequent n2 entries P2, last Hk entries Pk, H = 
[Hi,...,Hx] G C^^" and X = [X]", . . . ,xT ]T g c"><^. 
By convention, we shall assume Pi < . . . < Pk- 

Remark 1: The statement that \/iVH, X and W have 
independent entries of finite fourth order moment is meant 
to provide as loose assumptions as possible on the channel, 
signal and noise properties. In the simulations of Section V, the 
entries of H, W are taken Gaussian. Nonetheless, according 
to our assumptions, the entries of X need not be identically 
distributed, but may originate from a maximum of K dis- 
tinct distributions. This translates the realistic assumption that 
different data sources may use different symbol constellations 
(e.g., Af-QAM, M-PSK); the finite fourth moment assumption 
is obviously verified for finite constellations. 

Our objective is to infer the values of the powers 
Pi , . . . , Pk from the realization of the random matrix Y. This 
is the subject of Section IV. In the sequel, we introduce tools 
from large dimensional random matrix theory and we provide 
a thorough analysis of the eigenvalue distribution of -jgYY*^ 
as A'', n and M grow large at the same rate. 

III. Spectral analysis 

We start by analyzing the eigenvalue distribution of YY^ 
when n. A'' and M grow large at a similar rate. This is a 
fundamental prior step to the proper estimation of Pi , ... , Pk- 

A. Limiting spectrum of -^YY"^ 

We first define the Stieltjes transform of a (cumulative) 
distribution function. 

Definition 1: Let P be a distribution function. For z C\ 
M"*", the Stieltjes transform m{z) of F is defined as 



(4) 



For X G M a continuity point of F, we have the inverse 
Stieltjes transform formula 

1 



F{x) = — lim 

TT j/->-0+ 



•/ — < 



[m{t + iy)]dt. 



(5) 



In this section, we prove the following result 



Theorem 1: Let B 



N 



h-YY^, with Y defined as in 



M 



(3). Then, for M, N, n growing large with limit ratios 
M/N c, N/nk Cfe, < c,ci,...,ck < 00, the 
eigenvalue distribution function P^" of Bat, referred to as 
the empirical spectral distribution (e.s.d.) of Bat, converges 
almost surely to the deterministic distribution function P, 
referred to as the limit spectral distribution (l.s.d.) of B^f, 
whose Stieltjes transform mpiz) satisfies, for z G C"*", 



mF{z) = cmF{z) + (c — 1)- 

z 



(6) 



where mp_{z) is the unique solution with positive imaginary 
part of the implicit equation in m,p 



Pkf 



in which we denoted / the value 

/ = (1 — c)m,]7^ — czmpi 



(7) 



(8) 



The rest of this section is dedicated to the proof of Theorem 
1. First remark that (3) can be further simplified into 



HP5 



(9) 



Appending Y e C^^^ into the larger matrix Y e 



HP3 ah 




X 
W 



(10) 



we recognize that j^ YY ^ is a sample covariance matrix, for 
which the population covariance matrix ^ hph"+<7^I]v ^ is 
non-deterministic and the random matrix ( ^ ) has indepen- 
dent (non-necessarily identically distributed) entries with zero 
mean and variance 1. 

At this point, we need the following result. 

Proposition 1: Let Z„ S £^Ny.n j^^^g complex independent 
entries of zero mean, unit variance and finite 2 + e order 



moment, for some e > 0, and T„ G 



be Hermitian 



with e.s.d. converging almost surely to T, as A — > 00. Let 
A„ = ^Z„T„Zj^. Then, as n, A 00, N/n c > 0, the 
eigenvalue distribution of A„ converges weakly and almost 
surely to the distribution function A with Stieltjes transform 
mA{z), z G C+, being the unique solution with positive 
imaginary part of the equation in niA 

z = -— + -( -^—dT{t). (11) 
niA c J 1 + trriA 

Proof: The proof originates from Theorem 4.1 of [16] 
that states that, under the hypotheses of Proposition 1, the 
eigenvalue distribution of A„ converges weakly to some 
distribution function A whose Stieltjes transform mA{z) is a 
function of the Stieltjes transform of mT{z) and c only; mA{z) 
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is explicitly given by (4.4.4) of [16]. Now, in the special case 
where Z„ has independent and identically distributed (i.i.d.) 
entries of zero mean, unit variance and finite 2 + e moment, 
[10] and Theorem 4.3 of [16] show that mA{z) satisfies (11). 
But then, since mA{z) is only a function of c and T regardless 
of the distribution of the independent entries of Z„, mA{z) 
that solves (11) is the Stieltjes transform of A for the more 
general case. ■ 

Note that Proposition 1 can be equally stated when z e C~. 
In that case, mA{z) is the unique solution of (1 1) with negative 
imaginary part. 

From Proposition 1, since H has independent entries with 
finite fourth order moment, we have that the e.s.d. of HPH'^ 
converges weakly and almost surely to a limit distribution G 
as N,ni, . . . , riK — ^ oo with, N/n^ — Cfc > 0. For z E C+, 
the Stieltjes transform mciz) of G is the unique solution with 
positive imaginary part of the equation in tog. 



1 

TOG 



1 

^ Ck 1 

k=l 



Pk-niG 



(12) 



The almost sure convergence of the e.s.d. of HPH*^ en- 
sures the almost sure convergence of the e.s.d. of the matrix 
^hph"+o-2i„ 0^ Since tog(-z) evaluated at z e C+ is the 

Stieltjes transform of the l.s.d. of HPH'^ + ct^Iat evaluated 
at z + (7^, adding n zero eigenvalues, we finally have that the 
e.s.d. of ^ HPH"^+cr^iN ^ tends almost surely to a distribution 
H whose Stieltjes transform rnniz) satisfies 



mniz) = 



Co 



Co 



1 1 



1 + Co z ' 



(13) 



for z S C"*", where we denoted cq the limit of the ratio N/n, 
i.e., Co = (cj^^ + . . . + c^^)"^ 

As a consequence, the sample covariance matrix J- YY '^ 
has a population covariance matrix which is not deterministic 
but whose e.s.d. has an almost sure limit for increasing 
dimensions. Since X and W have entries with finite fourth 
order moment, we can again apply Proposition 1, and we have 
that the e.s.d. of = p-Y*^ Y converges almost surely to 
the limit F_ whose Stieltjes transform mp_{z) is the unique 
solution in C+ of the equation in mp_ 



1 1 / , 1 
+ - 1 + - 

mp C \ Cq 



t 



1 + tm F 



1 



1 + ^ 
cmp 



1 

nriF 



dH{t) 

1 

mp 



(14) 
(15) 



for all z e C+. 

For z E C+, raF_{z) E C+. Therefore — l/mi?(z) E C+ 
and one can evaluate (13) at — l/mi?(z). Combining (13) and 
(15), we then have 
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Fig. 2. Empirical and asymptotic eigenvalue distribution of -j^YY'^ when 
P has three distinct entries Pi = I, P2 = 3, P3 = 10, ni = n2 = na, 
CO = 10, c = 10, cr^ = 0.1. Empirical test: n = 60. 



where, according to (12), mG{^l/mp{z) — cr^) satisfies 
1 o 1 



mF_{z) 



tog(- 



mF_{z) 



^ 1 



Pk 



^ Cfc 1 + PfcTOG(- 



(17) 



m_F(z) 



Together with (16), this is exactly (7), with /(z) 



tog(- 



■) = (1 — c)toj7(z) — czmF{z) 



Since The eigenvalues of the matrices B n and B^ only 
differ by M~N zeros, we also have that the Stieltjes transform 
mp{z) of the l.s.d. of Bjv satisfies 



mp{z) — cmp{z) + (c ~ 1) — . 

~ z 



(18) 



This completes the proof of Theorem 1 . For further usage, 
notice here that (18) provides a simplified expression for 
tog(— 1/to_p(z) — ct^). Indeed we have. 



tog(— l/TOf^(z) — cr^) = —zmp{z)mp{z). 



(19) 



Therefore, the support of the (almost sure) l.s.d. F of Bjv 
can be evaluated as follows: for any z E C+, mp{z) is given 
by (6), in which mp_{z) is solution of (7); the inverse Stieltjes 
transform formula (5) allows then to evaluate F from mp{z), 
for values of z spanning over the set {z = x + iy,x > 0} 
and y small. This is depicted in Figure 2, when P has three 
distinct values Pi — 1, P2 = 3, P3 = 10 and ni = n2 — n^, 
N/n = 10, M/N = 10, ^ o.l, as well as in Figure 3 for 
the same setup but with P3 = 5. 

Two remarks on Figures 2 and 3 are of fundamental im- 
portance to the following. First, it appears that the asymptotic 
l.s.d. P of Bjv is compactly supported and divided into up 
to _ftr + 1 disjoint compact intervals, which we further refer 
to as clusters. Each cluster can be mapped onto one or many 
values in the set {a^, Pi, ... , Pk}- For instance, in Figure 3, 
the first cluster is mapped to cr^, the second cluster to Pi and 
the third cluster to the set {P2, Ps}- Depending on the ratios c 



5 



0.1 



0.075 



0.05 



0.025 



- Asymptotic spectiTjm 
Empirical eigenvalues 



H miim 




0.1 



3 5 
Estimated powers 



Fig. 3. Empirical and asymptotic eigenvalue distribution of -j^YY'^ when 
P has three distinct entries Pi = 1, P2 = 3, P3 = 5, ni = 712 = 113, 
CO = 10, c = 10, = 0.1. Empirical test: n = 60. 



and Co and on the particular values taken by Pi, ... , Pk and 
cr^, these clusters are either thin disjoint compact intervals, as 
in Figure 2, or they may overlap to generate larger compact 
intervals, as in Figure 3. We shall see, as is in fact required 
by the law of large numbers, that for increasing c and cq, 
the asymptotic spectrum tends to be divided into thinner and 
thinner clusters. The inference technique proposed hereafter 
relies on the separability of the clusters associated to each Pi 
and to (7^. Precisely, to be able to derive a consistent estimate 
of the transmitted power Pk, the cluster associated to Pk in F, 
number it cluster kp, must be distinct from the neighboring 
clusters {k — l)p and [k + \)f, associated to Pk-i and Pk+i 
respectively (when they exist), and also distinct from cluster 
1 in associated to . As such, in the scenario of Figure 
3, our method will be able to provide a consistent estimate 
for Pi, but (so far) will not succeed in providing a consistent 
estimate for either P2 or P3, since 2p = ?,p. We shall see that 
a consistent estimate for (P2 + P3)/2 is accessible though. 
Secondly, notice that the empirical eigenvalues of B^v are 
all inside the asymptotic clusters and, most importantly, in 
the case where cluster kp is distinct from 1, (fc — 1)_f and 
(fc + 1) p, observe that the number of eigenvalues in cluster kp 
is exactly Uk- This fact is referred to as exact separation. The 
exact separation for the current model originates from a direct 
application of the exact separation for the sample covariance 
matrix proven in [17] and is provided here in Theorem 3. This 
is further discussed in the subsequent sections. 



achievable if the cluster mapped to Pk in F is disjoint from 
all other clusters. The purpose of the present section is to 
provide sufficient conditions for cluster separability. 

To ensure that cluster kp (associated to Pk in F) is distinct 
from cluster 1 (associated to a^) and clusters ip, i 7^ fc, 
(associated to all other Pi), we assume now and for the rest 
of this article that the following conditions are fulfilled: (i) k 
satisfies Assumption 1, given as follows 

Assumption 1: 



^ (1 + P.ma.kf 



K 



E-- 



1 {PrrriGM+lY 



(1 + Prma^k+iY 



< 1 



(20) 
(21) 



with mG,i, ■ 
ma. 



, iTiG,K the K real solutions to the equation in 



^ 1 {Prmcf 



E-- 

tl^r (1 



(1 + Prmcy 



= 1 



(22) 



with the convention tuq = 0, 
and (ii) k satisfies Assumption 2 as follows. 
Assumption 2: Denoting, for j E {1, . . . , K}, 

Jg — ^ {i ^ j I * satisfies Assumption 1} , (23) 



I-cq {a'^mp^kGY 
Co (l + (r'^mp^ka)^ 



fr{ Cr{l + {x-^.^ + a^)mp^kaY 



r=kG + (^G,r + '^^)mF,kGy 

1- Co {a'^mp^ka+iY 
Co (1 + (j'^mp^ka+iy 



< c 



(24) 



E- 



(4,r+g')'^i,fco + l 

Kg 



(1 + (^G,r + CT2)TOj:,fcG + l) 

where Xq ^, Xq ■, i E {I, . . . , Kg}, are defined by 



< c 



(25) 



B. Condition for separability 

In the following, we are interested in estimating consistently 
the power Pk for a given fixed k £ {1,...,/^}. We recall 
that consistency means here that, as all system dimensions 
grow large with finite asymptotic ratios, the difference Pk — Pk 
between the estimate Pk of Pk and Pk itself converges to 
zero with probability one. As previously mentioned, we show 
by construction in Section IV that such an estimate is only 



with m 



G,l' "'ci' ■ 



1 



1 



'G,i 



K 

E 

r=l 
K 



1 



P, 



Cr 1 + Pr-m, 



G,i 



E-- 



p. 



P,.m, 



G,i 



(26) 



(27) 



(20), and mp^j, j E {1 



G,K, 



'■G Kg '^^G real roots of 
Kg + 1}, the j-th real root (in 
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increasing order) of the equation in mp 



100 



1 — Co {(j^mF_Y 



1 {x'^^. + a^yrriF; 



Co (1+0-2,7^,^)3 ^ (-j^ (-2.+ ^ ^ ^2),^^-)3 

1 (a^G^r + cr^)^™j^3 



+ 



Kg 



(28) 



Although difficult to fathom at this point of the article, the 
above assumptions will be clarified in the subsequent sections. 
We give hereafter a short intuitive explanation of the role of 
every condition. 

Assumption 1 is a necessary and sufficient condition for 
cluster kc, that we define as the cluster associated to Pf. in G 
(the l.s.d. of HPH'^), to be distinct from the clusters (fc— 1)g 
and {k+l)G, associated to Pk-i and Pk+i in G, respectively. 
Note that we implicitly assume a unique mapping between the 
Pi and clusters in G; this statement will be made more rigorous 
in subsequent sections. Assumption 1 only deals with the inner 
HPH'^ covariance matrix properties and ensures specifically 
that the powers to be estimated differ sufficiently from one 
another for our method to be able to resolve them. Note that, 
if Pi, ... , Pk are scaled by a common constant, then the 
solutions of (22) are scaled by the inverse of this constant; the 
separability condition is then a function of P2/^'i , • • ■ , Pk/ Pi 
and of the ratios ci , . . . , c/<- only. In Figure 4, we depict the 
critical ratio cq above which Assumption 1 is satisfied for all 
k, when K — 2 and ci = C2, as a function of P1/P2, i.e., the 
critical ratio cq above which the two clusters associated to Pi 
and P2 in G are disjoint. Observe that, as Pi gets close to P2, 
Co increases fast; therefore, to be able to separate power values 
with ratio close to one, an extremely large number of sensors 
is required. In Figure 5, the case A' = 3 is considered with 
ci = C2 = C3, Co — 10, and we let P2/P1 and P3/P1 vary; this 
situation corresponds to the scenarios previously depicted in 
Figures 2 and 3. Note that the triplet {Pi,P2,P3) = (1,3,5) 
is slightly outside the region that satisfies Assumption 1, and 
then, for this co, not all the clusters of G (and therefore 
of F) are disjoint, as confirmed by Figure 3. As for the 
triplet (1,3,10), it clearly lies inside the region that satisfies 
Assumption 1, which is sufficient to ensure the separability 
of the clusters in G, but not enough though to ensure the 
separability of the clusters in F. 

Assumption 2 deals with the complete Bat matrix model. 
It is however a non-necessary but sufficient condition so that 
cluster kp, associated to Pk in F, be distinct from clusters 
{k—l)p, {k+l)p and 1 (cluster 1 being associated to a^). The 
exact necessary and sufficient condition will be stated further 
in the next sections; however, the latter is not exploitable as 
is, and Assumption 2 will be shown to be an appropriate 
substitute. Assumption 2 is concerned with the value of c 
necessary to avoid (i) cluster ka (associated to Pk in G) to 
further spread on the clusters kc — 1 and kc + I associated 
to Pk-i and Pk+i and, more importantly, to avoid (ii) cluster 
1 associated to cr^ jjj p iq merge with cluster kp. As shall 
become evident in the next sections, when cr^ large, the 
tendency is for the cluster associated to to become large 
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Fig. 4. Limiting ratio cq to ensure separability of (Pi,P2). Pi < P2, 
K = 2,ci= C2. 



and spread over the clusters associated to Pi, then P2 etc. 
To counter this effect, one must increase c, i.e., take more 
signal samples. Figure 6 depicts the critical ratio c that satisfies 
Assumption 2 as a function of ct^, in the case K — 3, 
(Pi, P2, P3) - (1, 3, 10), Co = 10, ci = C2 = C3. Notice that, 
in the case c = 10, below ~ 1, it is possible to separate all 
clusters, which is compliant with Figure 2 where ct^ — g.l. 

As a consequence, under the assumption (proved later) 
that our proposed method cannot perform consistent power 
estimation when the cluster separability conditions are not met, 
we have two first conclusions: 

• if one desires to increase the sensitivity of the estimator, 
i.e., to be able to separate two sources of close transmit 
powers, one needs to increase the number of sensors (by 
increasing co), 

• if one desires to detect and reliably estimate power 
sources in a noisy environment, one needs to increase 
the number of sensed samples (by increasing c). 

In the subsequent section, we study the properties of the 
asymptotic spectrum of HPH'^ and in more details. 
These properties will lead to an explanation for Assumptions 
1 and 2. Under those assumptions, we shall then derive our 
novel power estimator 

IV. Multi-source power inference 

In this section, we prove our main result. 

Theorem 2: Let Bjy G C^^^ be defined as in Theorem 1, 
and A = (Ai, . . . , Xn), Ai < . . . < Aat, be the vector of the 
ordered eigenvalues of B^r. Further assume that the limiting 
ratios Cq, ci, . . . ,c/<-, c and P are such that Assumptions 1 
and 2 are fulfilled for some k £ {1, . . . , K}. Then, as N, n, 
M grow large, we have 



Pk - Pk 



where the estimate Pk is given by 







(29) 
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Fig. 6. Limiting ratio c as a function of to ensure consistent estimation 
of Pi = 1, P2 = 3 and P3 = 10, co = 10, ci = C2 = C3. 



if M ^ N, 



Pk = 



NM 



UkiM-N) 



(30) 



if M = N, 



N 



N 



UkiN-n) 



E E 



(31) 



in which Jiu = {Y.'l^i + ■ • ■ , Y.'l^i ^^i}' {m, ■ ■ ■ , Vn) are 
the ordered eigenvalues of the matrix diag(A) — i-y/A'v/A 
and (/ii, . . . , /ijv) are the ordered eigenvalues of the matrix 
diag(A) - ^%/A\/A^. 

Remark 2: We immediately notice that, if < n, the 
powers Pi, . . . ,Pi, with / the largest integer such that N — 
< 0' cannot be estimated. 



The approach pursued to prove Theorem 2 relies strongly on 
the original idea of [15]. From Cauchy's integration formula 
[18], 



Pk 



1 

2Td 



Pk-cj 



-duj 



1 



K 



Ck (p duj 

27ri ,/e Cr Pr — oj 



(32) 



for any negatively oriented contour 6^ C C, such that Pk 
is contained in the surface described by the contour, while 
for every i ^ k. Pi is outside this surface. The strategy is 
then the following: we first propose a convenient integration 
contour 6^ which is parametrized by a function of the Stieltjes 
transform mpiz) of the l.s.d. of B^v- We proceed to a variable 
change in (32) to express Pk as a function of niFiz)- We 
then evaluate the complex integral resulting from replacing the 
limiting mF{z) in (32) by its empirical counterpart mi?(z) = 

tr(Bjv — zIn)^^. This new integral, whose value we name 
Pk, is shown to be almost surely equal to Pk in the large N 
limit. It then suffices to evaluate Pk, which is just a matter of 
residue calculus [18]. 

We start by determining the integration contour Cfe. For 
this, we first need to study the distributions G and F in more 
details. 



A. Properties of G and F 

Let us introduce the following result on the l.s.d. of sample 
CO variance matrices, borrowed from [19] 

Proposition 2: Let A„ be defined as in Proposition 1 . Then 
the almost sure limiting Stieltjes transform m^(z) of the e.s.d. 
of A„, z € C+ admits a limit m\{x) when z a; G M*. If 
X is inside the support of A, then {x) is the only solution 
with positive imaginary part of the equation XA{ni) — x m 
the variable m, with XAim) defined, for — 1/m outside the 
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Fig. 8. XG{TnG) for mc real, P diagonal composed of three evenly 
weighted masses in 1, 3 and 5. Local extrema are marked in circles, inflexion 
points are marked in squares. 
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Fig. 9. xp{mp) for mp real, = 0.1, c = cq = 10, P diagonal 
composed of three evenly weighted masses in 1, 3 and 10. The support of F 
is read on the vertical axis. 



support of T, as 

11/" 

XA[m) = 1- - / T" 

m c J 1 



tm 



dT{t), 



(33) 



while, if X is outside the support of A, m°j^{x) is the only 
solution m of XAim) — x such that 2:^(771) > 0. Moreover, 
if for some m G M such that — 1/m is outside the support of 
T, x^(m) > 0, then XA{rn) is outside the support of A. 

The immediate corollary of Proposition 2 is that the com- 
plementary of the support Supp(yl) of A is the set {xAim)} 
for — 1/m outside the support of T such that a;^(m) > 0, 

Supp(^) = M \ {x I 3m e M, x = xa(to), x^(m) > 0} . 

(34) 

1} Support of G: First consider the matrix HPH'^, and 
let the function XGimc) be defined, for scalars mc e M* \ 



{-1/P,,...,-1/Pk}, by 



K 



XGimc) = 



niG 



r— 1 



PrTTlG 



(35) 



The function XQirao) is depicted in Figures 7 and 8, for 
the cases where cq — 10, c\~C2= C3 and (Pi, P2, -P3) equal 
respectively (1, 3, 10) and (1, 3, 5). As expected by Proposition 
2, xg(™g) is increasing for such that XGijno) is outside 
the support of G. Note now that the function xg presents 
asymptotes in the positions — 1/Pi, . . . , —1/Pk, 



lim XG(mG) 
lim XGirriG) 

'mG■n-^/P^) 



'OO, 



(36) 
(37) 



and that XGimG) — >■ 0+ as tug — > —00. Note also that, on 
its restriction to the set where it is non-decreasing, xg is 
increasing.^ To prove this, let and be two distinct 
points such that XG{niG) > and XG{rn*Q) > 0, and 
ttlq < ruG < 0, we indeed have,^ 



XGirriG) - XGirriG) 



■niG 



mGrriG 
K 



P2 



^ Cr iPr + —){Pr + ^) 



(38) 



Noticing that, for Pi > 0, 




Since we also have 



1 



^ 1 



p2 



1 



E 



1 



G ' 
P3 



> 



(42) 



> 0, (43) 



we conclude that the term in brackets in (38) is positive and 
then that xg{itig) — XG{rriG) > 0. Hence xg is increasing on 
its restriction to the set where it is non-decreasing. 

Notice also that xg, both in Figures 7 and 8, has exactly 
one inflexion point on each open set (— 1/Pi_i, — l/P^), for 

^we say here that a function f{x) is increasing if x < x* =^ fi^) ~ 
f{x*) > 0; if X < X* ^ /(^) ~ /(^*) ^ 0, we say that f{x) is non- 
decreasing. 

'this proof is borrowed from the proof of [15], with different notations. 
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i G {1, . . . , K}, with convention Pq = 0+. This is proven by 
noticing that XQirnc) = is equivalent to 



K 



1 



1 = 0. 



^ (1 + Prmaf 
Now, the left-hand side of (44) has derivative along mc. 



(44) 



p3 



{l + PrmaY' 



(45) 



which is always positive. Notice that the left-hand side of (44) 
has asymptotes for ma = —'i-/Pi for all i G {1, . . . , K}, and 
has limits as mc and 1/co — 1 as mc — oo. If 
Co > 1, Equation (44) (and then XQ{mG) = 0) therefore has a 
unique solution in (— 1/Pi_i, — for alH e {1, . . . , -ft"}. 
When xq is increasing somewhere on (— 1/Pi_i, — l/P^), 
the inflexion point, i.e., the solution to XQ{mG) = 0, in 
(— 1/P,_i, — 1/Pj) is necessarily found in the region where 
XG increases. If cq < 1, the leftmost inflexion point may not 
exist. 

From the discussion above and Proposition 2, it is clear that 

the support of G is divided into Kg < K compact subsets 
[xq j, Xq J, z e {1, . . . , Kg}. Also, if cq > 1, G has an addi- 
tional mass in of probability G{0) — G(0— ) = (cq — l)/co; 
this mass will not be counted as a cluster in G. Observe that 
every Pi can be uniquely mapped to a corresponding subset 
[xq j , Xq in the following fashion. The power Pi is mapped 
onto the first cluster in G; we then have 1g = 1- Then the 
power Pj is either mapped onto the second cluster in G if xg 
increases in the subset (— 1/Pi,— I/P2), which is equivalent 
to saying that a;^ (771,3 2) > for tog^2 the only solution 
to XQ{mG) = in (— 1/Pi, — I/P2); in this case, we have 
2(5 = 2 and the clusters associated to Pi and P2 in G are 
distinct. Otherwise, if xj^ (777,0,2) < 0, P2 is mapped onto the 
first cluster in F; in this case, 2g = 1. The latter scenario 
visually corresponds to the case when Pi and P2 engender 
"overlapping clusters". More generally, P, , j <E {1, . . . , K}, 
is uniquely mapped onto the cluster Jg such that 

JG = #{« < j I niin[xQ(mG^i),XQ(r77G,i+i)] > 0} , (46) 

with convention mG,K+i = 0, which is exactly 

JG = #{«<.? I i satisfies Assumption 1} , (47) 

when Co > 1. If co < 1, rnc.i, the zero of Xq in (— 00, — 1/Pi) 
may not exist. If co < 1, we claim that Pi cannot be evaluated 
(as was already observed in Remark 2). The special case when 
Co = 1 would require a restatement of Assumption 1 to handle 
the special case of Pi; this will however not be done, as it 
will turn out that Assumption 2 is violated for Pi if cr^ > 0, 
which we assume. 

In the particular case of the power P/. of interest in Theorem 
2, because of Assumption 1, XQ{mG,k) > 0. Therefore the 
index fee of the cluster associated to Pk in G satisfies kG = 
(fc— 1)g+ 1 (with convention Og = 0). Also, from Assumption 
1, x'Q{mG.k+i) > 0. Therefore {k + 1)g = fc^ + 1. In that 
case, we have that P^ is the only power mapped to cluster 
kc in G, and then we have the required cluster separability 
condition. 



2) Support of F: We now proceed to the study of F, the 
almost sure limit spectrum distribution of Bjv. In the same 
way as in the previous section, we have that the support of F 
is fuUy determined by the function X£{mp), defined for 
real, such that — 1 /mp_ lies outside the support of H, by 



XFimp) 



1 



+ 



1 



t 



CCq 



1 + tniF 



dH{t). (48) 



Figure 9 depicts the function xp_ in the case of Figure 2, 
i.e., X = 3, Pi = 1, P2 = 3, P;, = 10, ci = C2 = C3, Co = 10, 
c = 10, cr^ = 0.1. Figure 9 has the pecuhar behaviour that it 
does not have asymptotes as in Figure 7 where the population 
eigenvalue distribution was discrete. As a consequence, our 
previous derivations cannot be straightforwardly adapted to 
derive the spectrum separability condition. If co > 1, note also, 
although it is not appearing in the abscissa range of Figure 
9, that there exist asymptotes in the position mF_ = — l/a^. 
This is due to the fact that G(0) - G(O-) > 0, and therefore 
H (a^) — H [a^ —) > 0. We assume co > 1 until further notice. 

From Proposition 2, the support of F_ is complementary to 
the set of real normegative x such that x = xpiniF) and 
x'p{mp) > for a certain real mp, with x'p{mp) given by 

x'p{mF) = -^- / ,, dH(t). (49) 

— mF_ cco J (1 -h trriFj 

Reminding that H{t) = ■^p^G{t - ct^) -|- Y^5{t), this can 
be rewritten 

It is still true that xpirriF), restricted to the set of mp- where 
x'p{m£) > 0, is increasing. As a consequence, it is still true 
also that each cluster of H can be mapped to a unique cluster 
in F_. It is then possible to iteratively map the power Pk onto 
cluster kG in G, as previously described, and to further map 
cluster kG in G (which is also cluster fee in H) onto a unique 
cluster kp in F (or equivalently in P). 

Therefore, a necessary and sufficient condition for the 
separability of the cluster associated to Pfe in P reads 

Assumption 3: There exist two distinct real values 

(I) (r) 
''^F kr kG ^^^^ that 

1) ■'4Kk)>0,2;i.(rng,J>0 



2) there exist m^G k^''^%\ ^ ^ ^^^^ ^^'^^ ^ci'fn^Gk) ~ 

-l/m''2kG " ^cim'^GM) = -l/"^F,feG ~ 

satisfy 

b) and 

1 „ 1 



Pk-i < — 



<Pk< — 



<Pk+i (51) 



(0 ^ ^ ^ (r) 

^c'k ^c'k 
with the convention Pq = 0-|-, P;<-+i = 00. 
Assumption 3 states (i) that cluster kG in G is distinct from 
clusters {k — 1)g and (fc + 1)^ (Item 2b); this is another way 
of stating Assumption 1, and (ii) that the points m^pk^ — 
-l/(a;G(mg;,J + a2) and mg,^ 4 + 
(which lie on either side of cluster fcc in H) have respective 
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images x).'^ = XF{my,.J and xlj, = XFimYf^J by xp. 



such that Xp{mp\ ) > and x'p{mp',, ) > 0, i.e 



(r) 



xj^J, lie outside the support of F, on either side of cluster kp- 
However, Assumption 3, be it a necessary and sufficient 
condition for the separability of cluster kp, is difficult to 
exploit in practice. Indeed, it is not satisfactory to require 

(I) (r) 

the verification of the existence of such m.p and rrip j^^. 
More importantly, the computation of xf_ requires to know 
H, which is only fully accessible through the non-convenient 
inverse Stieltjes transform formula 



■F,kG> 



and 



H(x) = - lim 
TT y->-0 



f 

J — oc 



"m-H^t + iy)dt. 



(52) 



Instead of Assumption 3, we derive here a sufficient condi- 
tion for cluster separability in F_. Notice from the clustering 
of G into Kg clusters plus a mass at zero that (50) becomes 



x'pimp) = 



1 



mpf 
Co - 1 



1 Kg 



(1 + tmp)' 



cco (1 -|- a'^mp)'^ ' 



(53) 



where we remind that [xq j, J is the support of cluster i in 
G, i.e., a:;^ 1,2;;!; p . . • , a^G.Xc ' ^G.ifc images by xg 

of the 2Kg real solutions to Xg(mG) = 0. 

Observe now that the function — + tmp)'^, found in 
the integrals of (53), has derivative along t 



(1 + tmpy 



2t 



(1 + tmp)^ 



(1 + tmp) 



(54) 



and is therefore strictly increasing when < — 1/^ and 
strictly decreasing when > — For € (— 1/ (x^ ■ + 
cr^), —1/xq + a^), we then have the inequaUty 



x'p{mF) > 

Kg 

+ E 

r=i+l 



rap 




G,r 



+ cr- 



i(l + (4,, + Om^)2 



2\2 



(1 -I- {xq ,^ + u'^)mFy- 



+ 



Co - 1 



Co (1 + a^mp)'^ J ' 
(55) 



Denote fi{mp) the right-hand side of (55). Through the 
inequality (55), we then fall back on a finite sum expression 
as in the previous study of the support of G. In that case, 
we can exhibit a sufficient condition to ensure the separability 
of cluster kp from the neighboring clusters. Specifically, we 
only need to verify that fkG-i{''T^F,kG) > 0' with mp^kG the 
single solution to f'k^-iimp) = in the set (— kG-i~^ 

<^^)>-l/(2;G,fcG + ^^))' fkG{'mF,kG+i) > 0, with 

mp^kG+i unique solution to f'k {mp) = in the set 

(- V(3^G,feG +'^')' - V(%,fcG+i ))• ™^ ^'^^'^^y w^^t 

Assumption 2 states. 

Remember now that we assumed in this section cq > 1. 
If Cq < 1, then is in the support of H and therefore the 
leftmost cluster in F, i.e., that attached to a^, is necessarily 
merged with that of Pi. This already discards the possibility 
of spectrum separation for Pi and therefore Pi cannot be 



estimated. It is therefore not necessary to update Assumption 
1 for the particular case of Pi, when cq = 1. 

Therefore, Assumptions 1 and 2 ensure that (fc — 1)f < 
kp < {k+ l)p, kp ^ 1, and there exists a constructive way 
to derive the mapping k ^ kp. We are now in position to 
determine the contour 6^. 



B. Determination of Cfe 



From Assumption 2 and Proposition 2, there exist x^l^^ and 
outside the support of F, on either side of cluster kp, 
such that mpiyz) has hmits m^p = T^Fiix^kl) ''^f\g ~ 
mF°{x^i^^), as z — )■ x^iP^ and z x^^^^, respectively, with nip^ 
the analytic extension of mp in the points G 
M. These limits frSp and rn^p \^ are on either side of cluster 
fee in the support of ~l/H, and therefore — l/w^'^.^ — 

and —i/fn^pk^ — cr^ are on either side of cluster fee in the 
support of G. 

Consider any continuously differentiable complex path Tpk 

(I) (r) ' 

with endpoints and x)^^, and interior points of positive 
imaginary part. We define the contour Ci^^ as the union of 

''kf 



. and a;!'"-' € 



Tpk oriented from 



to xl 



i,f, and its complex conjugate 

T*p ^ oriented backwards from x).^ to x).^ . The contour C^,;; is 
clearly continuous and piecewise continuously differentiable. 
Also, the support of cluster kp in F_ is completely inside 
CF,fe, while the supports of the neighboring clusters are away 
from CiT^fc. The support of cluster fcc in H is then inside 
— llrnp{Qpk)f' and therefore the support of cluster kc in G is 
inside Cc.fe — —l/mF_{Gp^k) — f ^. Since mp_ is continuously 
differentiable on C \ R (it is in fact holomorphic there [19]) 
and has limits in .x^'^ and a;^'^\ ^G.k is also continuous and 
piecewise continuously differentiable. Going one last step in 
this process, we finally have that Pk is inside the contour 
C/c = — 1/7710(60, fc), while Pi, for all i 7^ fc, is outside 
Cfe. Since niG is also holomorphic on C \ M and has limits 



l/mF°{xl^^) - and -l/mF°{x 



a^, Qk is a 



in 

continuous and piecewise continuously differentiable complex 
path, which is sufficient to perform complex integration [18]. 

The contours Ci, 62, C3 originating from circular integration 
contours Gp^k of diameter [xfljx]^], k e {1,2,3}, for the 
case of Figure 2, are depicted in Figure 10. The points 



and x'"i^^ for kp € {1, 2, 3} are taken to be x)^,'^ — XF_{mip,kG), 
x^^^ = XFi{mF^kG+i)' with rriF^i the real root of /-{mp) = 
m\-l/{x+^^^+a^),-l/{xG .+a^)) when i € {1,2,3}, 
and we take the convention tuga = —1/(15 -|- cr^). 
Recall now that Pk was defined as 



(0 



Pk = Ck 



^ 1 



CO 



27ri ./Cfc ^ Cr Pr 



-du). 



(56) 



"^we slightly abuse notations here and should instead say that the support 
of cluster ko in H is inside the contour described by the image by — l/mj? 
of the restriction to C+ and C~ of Cf.A;, continuously extended to R in the 



points and -l/m£';.^. 



(r) 



11 




and 



13 10 

Real part of Cfc, fc = 1,2,3 

Fig. 10. Integration contours Cpi, Sf2 and Cf 3, for c = 10, cq = 10, 
Pi = 1,P2 = 3, Pa = 10. 

With the variable change a; = —l/mQ{t), this becomes 

K 



1 



mcit) 



E-- 

Cr 1 



Pr 



^ Cr 1 + Prmait) 
dt. (57) 



^ Co - 1 

Co / rnaitY 
From Equation (12), this simplifies into 

Pk = -^.i icotmcit) + (CO - 1)) ^^^^dt. (58) 

Using (16) and proceeding to the further change of variable 

t — ~l/mp_{z) — a^, (58) becomes 



/ 1 



\mp(z) 



zmF{z)mF{z) + 



Co - 1 



Co 



-mF_{z)mF{z) — zmF_{z)mp{z) — zmp_{z)m'p{z) 



Ck 
2711 



1 



co-1 



dz 



(l + a^mpiz)) + 
mi/{z) 



1 



Co zrnp{z) 



m'p{z) 



dz. 



(59) 



(60) 



zmpiz) mp{z)'^ mF{z)mp{z) 

This whole process of variable changes allowed us to 
describe Pk as a function of mp{z), the Stieltjes transform 
of the almost sure limiting spectral distribution of B^r, as 

— > oo. It then remains to exhibit a relation between Pk and 
the empirical spectral distribution of Bat for finite N . This is 
to what the subsequent section is dedicated to. 

C. Evaluation of Pk 

Let us now define rhp{z) and rhp[z) as the Stieltjes 
transforms of the empirical eigenvalue distributions of Bjv 
and B^, respectively, i.e., 



1 

^f{z) = — 



1 



(61) 



N M-N l 



(62) 



Instead of going further with (59), define Pk, the "empirical 
counterpart" of Pk, as 



n 1 



Pk - — : 



Uk Jqj^^ 



N 



1 



— (l + a^mpiz)) 
rn'p{z) 



N -n 



1 



n zmp^z) 



m'p{z) 



zmp{z) rhp{zY mp(z)mp{z) 



dz. 
(63) 



The integrand can then be expanded into nine terms, for 
which residue calculus [18] can easily be performed. De- 
note first rii,...,riN the N real roots of rhp{z) — and 
/ii, . . . , UN the N real roots of rhpiz) — 0. We identify three 
sets of possible poles for the nine aforementioned terms: (i) 
the set {Ai,...,AAr}n [x1l,x';^\ (ii) the set {771, ■ • • , ^7^} ^ 



,(0 J^) 



y^k!] (iii) the set {fii, . . . , (In} n [x]^'x\^},]. For 



M ^ N, the full calculus leads to 



Pk = 



NM 



nk{M~N) 



E 

l<i<N 
x<'' <rii<K!.'^' 



E 

l<i<N 

fee—"'— CI 



N 
nk 



N 
nk 



E 

l<i<Af 
E<'' <r]i<xi''^ 

fcc ' k IT 



E 

l<i<Af 

fcu" — — fcc- 



E 

l<i<N 



E 

l<i<Ar 
' <Ai<a:i 



(64) 



Details are given in Appendix A. Now, we know from 
Theorem 1 that rhp{z) mp[z) and rhp_{z) mF_{z) 
as — > 00. Observing that the integrand in (63) is uniformly 
bounded on the compact Qp,k, the dominated convergence 
theorem [20] ensures i\. Pk- 

To go further, we now need to determine which of 
Ai, . . . , \m, ?7i, . . . , ?7jv and ^1, . . . , /x^v lie inside Qp^k- This 
requires a result of eigenvalue exact separation that extends 
the earlier results of [21], [17], as follows 

Theorem 3: Let B„ = (l/n)T|x„X,^T| G C?'''?', where 
we assume the following conditions 

1) X„ e C^^" has entries Xij, l<i<p, l<j<n, ex- 
tracted from a doubly infinite array {xij} of independent 
variables, with zero mean and unit variance. 

2) There exist K and a random variable X with finite fourth 
order moment such that, for any x > 0, 



1 



J2 P{\x^J\> x) <KP{\X\> x) (65) 



for any ni, 112. 
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3) There is a positive function il){x) t oo as a; oo, and 
M > 0, such that 



maxE|x?-|^(|xij|) < M. 

ij 



(66) 



p the ordered eigenvalues of the 
^■"P. Then, we have that 



A) p = p{n) with c„ = p/n ^ c > as n ^ oo. 

5) For each n, T„ G C^^p is Hermitian nonnegative defi- 
nite, independent of {xij}, satisfying Hn = F'^'^ => H, 
H a nonrandom probability distribution function, almost 
surely. is any Hermitian square root of T„. 

6) The spectral norm ||T„|| of T„ is uniformly bounded 
in n almost surely. 

7) Let a,b > 0, nonrandom, be such that, with probability 
one, [a, b] lies in an open interval outside the support 
of F'^"'^" for all large n, with F^''^ defined to be the 
almost sure l.s.d. of (l/n)Xj^T„X„ when H = G and 
c = y. 

Denote > . . . > 
Hermitian matrix Y e 

1) P(no eigenvalue of B„ lies in [a, b] for all large n) = 
1. 

2) If c(l - H{0)) > 1, then .xo, the smallest value in the 
support of F"'^, is positive, and with probability one, 
A^" ^ a^o as n ^ oo. 

3) If c(l - H{0)) < 1, or c(l - (0)) > 1 but [a, b] is 
not contained in [0,xo], then mpa,H{a) < mpa,H{b) < 
0. With probability one, there exists, for all n large, 
an index > such that A^" > —l/mpc,H{b) and 
Xj":^_i > —l/mpc,H{a) and we have 

P{Xf^" > b and AJ^-^^ < a for all large n) = 1. (67) 

Theorem 3 is proven in Appendix B. 

To apply Theorem 3 to B^y in our scenario, we need to 
ensure all assumptions are met. Only Items 2-6 need particular 
attention. In our scenario, the matrix X„ of Theorem 3 is 
(^), while T„ is T ^ (^HPH"+^=^i„ 0^ The latter has 
been proven to have almost sure l.s.d. H, so that Item 5 is 
verified. Also, from the result of [21] upon which Theorem 3 is 
based, there exists a subset of probability one in the probability 
space that engenders the T over which, for n large enough, 
T has no eigenvalues in any closed set strictly outside the 
support of H; this ensures Item 6. Now, from construction, X 
and W have independent entries of zero mean, unit variance, 
fourth order moment and are composed of at most K + 1 
distinct distributions, irrespectively of M. Denote Xi,. . . , X^, 
d < K + 1, d random variables distributed as those distinct 
distributions. Letting X = |Xi| + . . . -|- \Xd\, we have that 

V Pi\zij\>x)<p(y2\Xi\>x] (68) 

= P{\X\ > x), (69) 

where Zij is the {i,jy*^ entry of ( ). Since all X^ have finite 
order four moments, so does X and Item 2 is verified. From 
the same argument, Item 3 follows with 4>{x) = x^. Theorem 
3 can then be apphed to B^y. 

The corollary of Theorem 3 applied to B;y is that, with 
probability one, for N sufficiently large, there will be no 



eigenvalue of Bjv (or B^y) outside the support of F, and 
the number of eigenvalues inside cluster kp is exactly rife. 
Since &p^k encloses cluster kp and is away from the other 
clusters, {Ai, . . . , Ajv} n [x^H^x^^^] ^ e Jik] almost 

surely, for all large N . Also, for any i e {1, . . . , N}, it is 
easy to see from (61) that rnp{z) — > oo when z '\ Xi and 
mp{z) — > — oo when z I Xi. Therefore 771^(2) = has at least 
one solution in each interval (Ai_i,Ai), with Ao = 0, hence 
^1 < Ai < ^2 < • • • < ^J'N < Xn- This imphes that, if fco is 
the index such that Cp^k contains exactly Xkg, ■ • . , AfcQ+(„^_i), 
then Gp^k also contains {fiko+i, ■ ■ ■ , ^J'ko+ink-^^)}■ The same 
result holds for rjkg+i, ■ ■ ■ ,r]k„+(nk~i)- When the indexes 
exist, due to cluster separability, rjkg-i and Hko-i belong, for 
N large, to cluster kp — 1. We are then left with determining 
whether Hk^ and rjko are asymptotically found inside 6^,^. 

For this, we use the same approach as in [15], by noticing 
that, since is not included in Cfe, one has 

— / -duj = 0. (70) 

Performing the same changes of variables as above, we have 
that 

—mF_(z)mp{z) — zmF_{z)mF{z) — zmF_{z)m'p{z) 



= 0. 



z'^mp_{zYmp{zY 



dz 



(71) 



For N large, the dominated convergence theorem ensures 
again that the left-hand side of the (71) is close to 

—mF_{z)rhp{z) — zrh'p{z)mp{z) — zrhp^{z)m'p{z) 



ef.,k z'^rhF_{zYrhp{zY 
Residue calculus of (72) then leads to 



E 



E 1- E 



l<i<N 



l<i<Af 



l<i<N 



Since the cardinalities of {i,rii e [Xkl t ^''j^p]} and {i,/tii € 



-dz. 
(72) 



0. 



(73) 



X 



(I) Jr) 



,x^^J^]} are at most nk, (73) is satisfied only if both 
cardinalities equal Uk in the limit. As a consequence, ^kn G 



[x^kl'4l] and Vko G [4^,Xkl]- For N large, iV ^ M, this 
allows us to simplify (64) into 

NM 



Pk 



nk{M-N) ^ 



E 

<i<N 



(74) 



N. 



with probability one. The same reasoning holds for M 
This is our final relation. 

It now remains to show that the ryj and the jii are the 
eigenvalues of diag(A) — ;^\/A\/A and diag(A) — -p-\/A\/A 
respectively. For this, we need the following lemma. 

Lemma 1: Let A G R^^^ be diagonal with entries 
Ai , . . . , Ajv, and let y e M^. Then the eigenvalues of A— yy^ 
are the N real solutions of the following equation in x. 



Vt 



N 

^'Xi-x 
1=1 



1. 



(75) 
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Proof: Let A be an eigenvalue of A — yy"^ . For a certain 
non-zero vector x G C^, we then have the equivalent relations 

(A-yy'^)x = Ax, (76) 

(A - AIiv)x = y^xy, (77) 

x = yHx(A-AI^)-V, (78) 

y"x = y'^xy'^(A-AlAr)-V, (79) 

l = y"(A-AIjv)-V- (80) 

Since A is diagonal, denoting ej e the vector such that 
^i,3 = <^i. we finally have 



^ A,; -A 



= 1. 



(81) 



i=l 



Applying Lemma 1 to A = diag A and y = y A, we find 
that the eigenvalues of diag(A) — -^^/X^/X are the solutions 



of 



N i_ 
N' 



E 

i=l 



A,; - X 



which is equivalent to 



t — 1 



0, 



(82) 



(83) 



whose solutions are by definition r]i,. . . ,r]N- The same argu- 
ment applies similarly to Hi,. . . , Hm- Incidentally, this remark 
was already noticed in [22]. 

We end this section by a short discussion on the conse- 
quences of Theorem 2. 

D. Discussion 

Theorem 2 states that, under spectrum separabihty condition 
for all Pfe, fc e {1, . . . , K}, when ni , . . . , uk are known a pri- 
ori to the receiver, then Pi, ... , Pk are consistent estimators 
for Pi, ... , Pk- Now, in practice, it is rare that rii, . . . , uk 
and even K are, a priori known to the receiver. However, if 
separability is assumed, then one can estimate simultaneously 
K,ni, . . . , uk and Pi, ... , Pk- This is performed by (i) 
determining the clusters of the empirical eigenvalues of Bjv, 
which determines K, (ii) counting the number of eigenvalues 
in each cluster to determine the multiphcities ni , . . . , uk and 
(iii) evaluating Pi , • • • , Pk from Theorem 2. 

However, step (i) may not be obvious. In particular, when 
the total number n of transmit antennas is small, when the 
typical cluster size is large or when the inter-cluster spacing 
is small, it is non-trivial to determine what eigenvalues form a 
cluster. To solve this critical issue, studies are being currently 
carried out that aim to determine second order statistics of 
xhanks to second order statistics on F^™, it will be 
possible to design estimators of Pi , . . . , Pk that take into 
account the probability of Bjv being an appropriate model 
for the estimated Pi , . . . , P^. for every hypothesis K for the 
number of transmit source and every hypothesis (ui, . . . ,nf.) 
for the number of antennas for each of these sources. We 
hereafter provide an alternative ad-hoc technique to partially 



solve the problem of determining K and ni , . . . , uk based on 
Theorems 1 and 2. 

In the following, we assume for readability that we know 
the number K of transmit sources (taken large enough to cover 
aU possible hypotheses), some having possibly transmit 
antermas. The approach consists in the following steps: 

1) we first identify a set of plausible hypotheses for 
ni, . . . , riK- This can be performed by inferring clusters 
based on the spacing between consecutive eigenvalues: 
if the distance between neighboring eigenvalues is more 
than a threshold, then we add an entry for a possible 
cluster separation in the list of all possible positions of 
cluster separation. From this list, we create all possible 
ii'-dimensional vectors of eigenvalue clusters. Obvi- 
ously, the choice of the threshold is critical to reduce 
the number of hypotheses to be tested; 

2) for each if-dimensional vector with number of antennas 
ni,...,nK, we use Theorem 2 in order to obtain 
estimates of the Pi, ... , Pk (some being possibly null); 

3) based on these estimates, we compare the e.s.d. P^" of 
Bat to the distribution function F defined as the l.s.d. of 
the matrix model Y = HPX -I- W with P the diagonal 
matrix composed of hi entries equal to Pi, "-2 entries 
equal to P2 etc. up to fiK entries equal to Pk- The 
comparison can be performed based on different metrics. 
In the simulations carried hereafter, we consider as a 
metric the mean absolute difference between the Stieltjes 
transform of P^" and of F on the segment [—1, —0.1]. 

Note that the above process can bring an interesting feature 
linked to the cluster separabihty problem discussed along 
this article. Indeed, if two subsequent powers P; and Pj+i 
are close to one another, then the separability condition of 
Assumptions 1 and 2 is not verified. If one knows rij and 
rij+i and blindly uses the estimator of Theorem 2, the result 
can be catastrophic as the estimator is unreliable. On the 
contrary, if m and rii+i are unknown and one uses the above 
process, it is very Ukely that the distinct sources with close 
power will be assumed to be a single source with power 
equal to the estimate of (Pj + Pi+i)/2 and embedded with 
rij -|- rij+i antennas. For practical blind detection purposes 
in cognitive radios, this leads the secondary network to infer 
a number of transmit entities that is less than the effective 
number of transmitters. In general, this would not have serious 
consequences on the decisions made by the secondary network 
but this might at least reduce the capabihties of the secondary 
network to optimally overlay the licensed spectrum. Further 
work is also being carried out to go past the cluster separability 
assumption; specifically, methods for estimating the number of 
Pi associated to any cluster jp are under study. 

V. Simulations 

In this section, we provide simulation results to assess 
the performance of Theorem 2 when K, and ni , . . . , uk 
are known, to compare this performance against alternative 
estimation methods and finally to evaluate the performance 
of the ad-hoc approach discussed in Section IV-D. In order 
to underUne some precise features of the advantages of our 
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novel method, we will use two simulation models. The first 
model, already presented in Figure 2, involves a scenario with 
clear separation between clusters, while the second model will 
consider the case of co-located clusters. 

The estimator of Theorem 2 wiU be compared against two 
methods, which we describe below. 

A. Alternative methods 

1 ) Strongly consistent estimator for M ^ N and N ^ n: 
The first method is the classical estimator that assumes that the 
sample dimension M is much larger than the sensor dimension 
N, while N is much larger than the source dimension n. In this 
case, it is easy to see that the e.s.d. of B^r tends to a mass in 
<j^. However, the first n eigenvalues of Ba? are asymptotically 
greater than cr^ and it is also clear that the e.s.d. of the 
projection of Bat on the eigenspace associated to its largest n 
eigenvalues tends to K masses in Pi + a^, . . . , Pk + cr^. This 
leads to the strongly consistent estimator of Pk given by 



with 



^ N-n 



(84) 



i=l 

and we recall that Ai < . . . < Ajv are the eigenvalues of Bjv. 
The strong consistence is with respect to the rates n — >■ oo, 
N/n 00 and M/N oo. Note that we take an estimator 
for cr^ instead of cr^ itself in order to be coherent with Theorem 
2 which does not require any a priori information on ct^. We 
will refer to this estimator as the classical method. 

2) Estimator based on strongly consistent moment esti- 
mates: The second method is a technique issued from free 
probability theory, which is based on moments of the l.s.d. 
of B AT . As such, we will refer to this method as the moment 
method. It consists in computing the first moments of the e.s.d. 
of Bat, i.e., i tr {j^YY^f, for fc = 1, . . . , if, from which 
the deconvolved moments ^(niPf + . . . -|-nifP^) of can 
be evaluated, see e.g., [23]. These estimated moments can be 
expressed as polynomials of the moments of F^", which is 
convenient from a practical point of view although it leads to 
serious shortcomings in terms of estimator accuracy. Indeed, 
small deviations in the low order moments of f around 
the corresponding moments of F lead to large deviations in 
the estimation of the high order moments of . 

One can then retrieve the vector . . . , 

whose distribution function has for first K moments the first K 
estimated moments of F^ . This is performed using Newton- 
Girard polynomial formulas [24], which boils down to finding 
the roots of a polynomial of order K. The value p^™°™) is 
the estimate of Pk- Computing p^™"™^ requires in particular 
that if , ni , . . . , TiK and are known. The main shortcoming 
of the Newton-Girard inversion is that the polynomial to be 
solved may have purely imaginary roots. This issue, added to 
the deviations in the estimated moments, contribute to rather 
poor estimation accuracies unless the system dimensions are 
very large. However, as opposed to the classical method and 



the novel Stieltjes transform approach, the moment method 
does not require any assumption of cluster separabihty to be 
valid. 



B. Results 

1 ) Cluster separability limit: We start with a demonstration 
of the performance of the novel estimator with respect to the 
satisfaction of the cluster separability assumption. We consider 
the model presented in Figure 2, i.e., K = 3, Pi = 1, P2 = 3, 
P3 = 10, ni/n = n2/n = m/n = 1/3 and n/N = N/M = 
1/10. The SNR, defined as SNR = l/cr^, ranges from -15 dB 
to 20 dB. The entries of X are QPSK-modulated and those of 
H and W are Gaussian distributed. In Figure 11, we present 
simulation results in terms of normahzed mean square error 
(NMSE) in the estimates of the individual Pk, both for n = 60 
and n = 6. For future need, we define this system model with 
n = 6 as Scenario (a). The NMSE for power Pk is given by 



NMSE = E 



p2 



(85) 



where the expectation is taken over the random reahzations of 
the matrices H, X and W. 

Note how steep the mean square error curves increase below 
a given SNR value. This intuitively corresponds to the tipping 
point where the cluster separability assumptions are no longer 
verified. Especially here, this corresponds to the point where 
Assumption 2 no longer holds. Now, remembering the results 
of Figure 6, observe that the horizontal line c = 10 crosses 
the respective curves of validity of Assumption 2 around the 
SNR values where Figure 11 shows steep curve increase. 
This indicates that our novel estimator is indeed inappropriate 
when Assumption 3 is not satisfied. This also validates the 
accuracy of Assumption 2, which we recall is only a sufficient 
condition for cluster separabihty. Note also that, as long as 
cluster separation is achieved, the performance of the Stieltjes 
transform algorithm goes quickly down to a constant level 
(with respect to the SNR) which is a function of the amphtude 
of the values of n, N and M. 

2) Performance comparison: We first compare the classical 
method against the novel Stieltjes transform approach for 
Scenario (a). Under the hypotheses of this scenario, the ratios 
c and Co equal 10, leading therefore the classical detector to 
be almost asymptotically unbiased. We therefore suspect that 
the NMSE performance for both detectors is alike. This is 
described in Figure 12, which suggests as predicted that in 
the high SNR regime (when cluster separability is reached) the 
classical estimator performs similar to the Stieltjes transform 
method. However, it appears that a 3 dB gain is achieved by the 
Stieltjes transform method around the position where cluster 
separability is no longer satisfied. This translates the fact that, 
when subsequent clusters tend to merge as increases, the 
Stieltjes transform method manages to track the position of the 
powers Pk while the classical method keeps assuming each 
Pk is located at the center of cluster hp- This observation is 
very similar to that made in [25], where an improved MUSIC 
estimator is introduced that pushes further the SNR position 
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Fig. 11. Normalized mean square error of individual powers Pi, P2, P3, 
Pi = 1,P2 = 3,P3 = 10, ni/n = n2/n = ns/n = 1/3 ,n/N = 
N/M = 1/10, for 10, 000 simulation runs. 
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Fig. 12. Normalized mean square error of individual powers Pi, P2, P3, 
Pi = 1,P2 = 3,P3 = 10, rii/n = n2/n = ns/n = 1/3 ,n/N = 
N/M = 1/10, 71 = 6. Comparison between classical and Stieltjes transform 
approach. 

where the performance of the classical MUSIC estimator 
decays significantly. 

We now consider another model, for which the classical 
estimator is largely biased. We now take K = 3, Pi — 1/16, 
P2 = 1/4, P3 = 1, ni/n — n2/n ~ n-s/n ~ 1/3 and 
n = 12, iV = 24 and M = 128. The entries of X are 
still QPSK-modulated while the entries of H and W are still 
independent standard Gaussian. This model is further referred 
to as Scenario (b). We first compare the performance of the 
classical, Stieltjes transform and moment estimators for an 
SNR of 20 dB. Figure 13 depicts the distribution function of 
the estimated powers in logarithmic scale. The Stieltjes trans- 
form method appears here to be very precise and seemingly 
unbiased. On the opposite, the classical method, with a slightly 
smaller variance shows a large bias as was anticipated. As for 



Fig. 13. Distribution function of the estimators P^, Pfe, ^nd 
for k G {1,2,3}, Pi = 1/16, P2 = 1/4, P3 = 1, ni = n2 = na = 4 
antennas per user, A'^ = 24 sensors, M = 128 samples and SNR = 20 dB. 
Optimum estimator shown in dashed lines. 



the moment method, it shows rather accurate performance for 
the stronger estimated power, but proves very inaccurate for 
smaller powers. This entails from the inherent shortcomings 
of the moment method. The performance of the estimator Pj, 
will be commented in Section V-B3. 

We then focus on the estimate for the larger power P3 and 
take now the SNR to range from —15 to 30 dB under the 
same conditions as previously and for the same estimators. 
The NMSE for the estimators of P3 is depicted in Figure 14. 
The curve marked with squares will be commented in Section 
V-B3. As already observed in Figure 13, in the high SNR 
regime, the Stieltjes transform estimator outperforms both 
alternative methods. We also notice the SNR gain achieved by 
the Stieltjes transform approach with respect to the classical 
method in the low SNR regime, as already observed in Figure 
12. However, it now turns out that in this low SNR regime, 
the moment method is gaining ground and outperforms both 
cluster-based methods. This is due to the cluster separability 
condition which is not a requirement for the moment approach. 
This indicates that much can be gained by the Stieltjes 
transform method in the low SNR regime if a more precise 
treatment of overlapping clusters is taken into account. 

3) Joint estimation of K, rik, Pk-' So far, we have assumed 
that the number of users K and the number of antennas 
per user were perfectly known. As discussed in Section 
IV-D, this may not be a strong assumption if it is known by 
advance how many antennas are systematically used by every 
source or if another mechanism, such as in [9], can provide 
this information. Nonetheless, these are in general strong 
assumptions to take. Based on the ad-hoc method described 
in Section IV-D, we therefore provide the performance of our 
novel Stieltjes transform method in the high SNR regime when 
only n is known; this assumption is less stringent as in the 
medium to high SNR regime, one can easily decide which 
eigenvalues of Bat belong to the cluster associated to cr^ 
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Fig. 14. Normalized mean square error of largest estimated power P3, 
Pi = 1/16, P2 = 1/4, P3 = 1, m = n2 = na = 4 ,Af = 24, M = 128. 
Comparison between classical, moment and Stieltjes transform approaches. 



and which eigenvalues do not. We denote the estimator 
of Pk when K and are unknown. We assume 

for this estimator that all possible combinations of 1 to 3 
clusters can be generated from the n — 6 observed eigenvalues 
in Scenario (a) and that all possible combinations of 1 to 
3 clusters with even cluster size can be generated from the 
n = 12 eigenvalues of B^v in Scenario (b). For Scenario 
(a), the NMSE performance of the estimators and is 
proposed in Figure 15 for the SNR ranging from 5 dB to 30 
dB. For Scenario (b), the distribution function of the inferred 
is depicted in Figure 13, while the NMSE performance for 
the inference of P3 is proposed in Figure 14; these are both 
compared against the classical, moment and Stieltjes transform 
estimator. We also indicate in Table V-B3 the percentage of 
correct estimation of the triplet (ni, n2, Jt-s) for both Scenario 

(a) and (b). In Scenario (a), this amounts to 12 such triplets 
that satisfy tij. > 0, rii + n2 + = 6, while in Scenario 

(b) , this corresponds to 16 triplets that satisfy rik G 2N, 
ni + 122 + = 12. Observe that the noise variance, assumed 
to be known a priori in this case, plays an important role 
with respect to the statistical inference of the n^. In Scenario 
(a), for an SNR greater than 15 dB, the correct hypothesis 
for the rik is almost always taken and the performance of 
the estimator is similar to that of the optimal estimator. In 
Scenario (b), the detection of the exact cluster separation is 
less accurate and the performance for the inference of P3 
saturates at high SNR to -16 dB of NMSE, against -19 
dB when the exact cluster separation is known. It therefore 
seems that in the high SNR regime the performance of the 
Stieltjes transform detector is loosely affected by the absence 
of knowledge about the cluster separation. This statement is 
also confirmed by the distribution function of P^ in Figure 13, 
which still outperforms the classical and moment methods. We 
underline again here that this is merely the result of an ad-hoc 
approach; this performance could be greatly improved if e.g., 
more is known about the second order statistics of P^" . 



SNR 


RCI (a) 


RCI (b) 


do 


U.O^ t 


u. looy 


10 dB 


0.9026 


0.4798 


15 dB 


0.9872 


0.4819 


20 dB 


0.9910 


0.5122 


25 dB 


0.9892 


0.5455 


30 dB 


0.9923 


0.5490 



TABLE I 

Rate of correct inference (RCI) of the triplet (ni, 712,713) for 

SCENARIOS (A) AND (B). 
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Fig. 15. Normalized mean square error of individual powers Pi, P2: ^"3 



,n/N = N/M 



Pi = 1, P2 
= 1/10, n = 



3, P3 = 10, ni/n = n^/n - 
, 10, 000 simulation runs. 



: 713/71 = 1/3 



VI. Conclusion 



In this paper, a blind multi-source power estimator was 
derived. Under the assumptions that the ratio between the 
number of sensors and the number of signals sources is not too 
small and the source transmit powers are sufficiently distinct 
from one another, we derived a method to infer the individual 
source powers if the number of sources is known, which 
was shown to outperform alternative estimation techniques in 
the medium to high SNR regime. We then briefly discussed 
the joint estimation of the number of transmit sources, the 
number of antennas of each source and the transmit powers, 
which appeared in simulation to perform well in the high 
SNR regime. The novel method is moreover computationally 
efficient and is particularly robust to small system dimensions. 
As such, it is particularly suited to the blind detection of 
primary mobile user in future cognitive radio networks. 
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Appendix A 
Residue calculus 

The integrand of Pk in (63) can be expanded as 

N 1 N rh'p{z) N m'p{z) 



n zmF_{z) n mF_{zY n mF^{z)inF{z) 
N-n 1 N-n rhpiz) 



n z'^rhF_{z)rhF{z) n zrhF_{zYrhF{z) 
N-n rh'p{z) Na^l Na^ m'p{z) 

n z n mF_{z) 



n z'mF{z)^rhF_{z) 
Na'^ rh'p{z) 



(86) 



n mpiz) 

First assume the case M ^ N. Numbering the nine terms in 
order, we have that (1) has poles in z e {ryi, . . . , ryjv}, where 
mF^z) = 0. Applying I'Hospital rule, all poles have order 1 
and the corresponding residues are 

N z-r]i N 1 



lini - 

z^rji n zmp(z) 



n riim'p{rii) 



(87) 



As for (2), it is the derivative of —l/7fiF_{z), which is well- 
behaved inside GF,k, so it does not have poles. The term (3) 
has poles of order 1 in 2; e {r/i, . . . , Jjjv} as well and we have 
the residue 



lim 



N {z - T]i)m'p{z) ^ N N m'pjijk) 
z->-vi n 7hFiz)mF{z) n M — N m'p{r}k) 

the last equaUty being obtained from the fact that 

^ , , M ^ , , M-Nl 
mF{z) = ^rriFiz) + 

and fhp{rii) = 0. It also has poles in z e {/^i, . . . , /xjv}, where 
rhF{z) = 0. These are order 1 poles and we have 



% (88) 



(89) 



N {z - Hi)m'p{z) N M 

1™ / ^ - / N = —17 — 

z^tii n mp[z)mp[z) n M — i\ 



(90) 



Term (4) is shown in a similar way to have residues 
-"^l^^^'M and ^i^^j-ra'pi^,,) , i e 
{!,..., AT}. Term (5) has residues ^^M^Vi and 

— {M-Nr m'^T^.) ' * ^ {1,...,N}. Term (6) has 



residues — 



n {M-N}'-' m'p{rii) 



and 



N-n MN 



{1, . . . , N}. Term (7) has a pole in z = but we already 
know that is not inside &F,k, so this is already discarded. 
Term (8) has poles in 2; e {Ai, . . . , Ajv} of residue and 
poles in z S {771, . . . ,?7Ar} of residue Similarly term 

(9) has poles in 2; e {Ai, . . . , Ajv} of residue and poles 
in 2 e {ni, . . . , hn} of residue 

Summing together the 9 terms and remarking that 



N 



-m'p{z) 



M 



M-N 



rh'p{z) 



1 



(91) 



M-N 
we obtain exactly (64). 

Assume now M = N, in which case 771^(2) = mF{z). 
It can be readily seen that the terms (4) to (6) are the 
derivative of — -^^^ ,2^^ so that they have residue 0. 

Th Z III p yZ J 

The only remaining term here is (1), whose residues are the 

jv 1 

n ■nim'p(rii) ' 



Appendix B 
Proof of exact separation 

Theorem 3 is a generalization from the assumption of iden- 
tical distribution of the Xij 's, the proof of which is contained in 
the two papers [21] and [17], which, with some modifications, 
appear as Chapter 6 in [16]. The proof uses previous articles 
that need to be updated as well. We shall therefore go through 
the necessary steps that need to be modified, taking reference 
to all papers successively. 

We shall assume for simplicity in the following that the 
matrices T„ are deterministic, converges in distribution to 
H and that ||T„|| is uniformly bounded. The generalization 
to random T„ follows from TonelU's theorem [20]. Indeed, 
let X be the probability space that engenders the X„ and 
T the probability space that engenders the T„. Let A be 
any of the events in Conclusion 1), 2), or 3), claimed to 
occur with probability one. Assume that Theorem 3 holds for 
deterministic T„ satisfying Assumptions 5), 6), and 7). Let 
t G T be an element of the intersection of these events. Then 
lA{t,x) = 1 for all X contained in a subset of X having 
probability one. Therefore, by Tonelli's theorem, denoting 
0" X X the product space of T and X, we have that 



/ 

-III 



lA{t,x)dP'jxxit,x) 



lA{t,x)dPx{x) 



dP7{t) = 1, 



(92) 



and Theorem 3 therefore holds true if it holds true for T„ 
deterministic. 

A. Extension of [26] 

The first step is to extend the work in [26] on the largest 
eigenvalue of S„ = iX„XjJ, where X„ = [xij) is p x n, 
p = p{n), and p/n — >^ y > as n ^ 00. Checking the 
assumptions in [26], we change the six conditions of Page 
518 to 



(1) X, 



1, 2, . . . , p; j = 1, 2, . . . , n are independent for 



each n, 

(2) I I < ?7„A/n, where ??„ i 0, 

(3) ^X^J = 0, 

(4) EI4I < 1, 

(5) E|xyf < {r^n^Y'\ for/> 2, 

(6) E|4|<c(?7„v^)'-3, for;>3. 

By the same argument given there (no difference for com- 
plex random variables), the inequality 

Etr(S„)'= < T?'^ (93) 

holds for any r] > h = [1 + ^Jy)"^ provided k is chosen such 
that 

(a) k/ logn — )• 00, 

(b) r/lfc/logn 0. 

This impHes that P(Amax(S„) > h + e) = o{n~*) for any 
given e > and t > 0. 

Remark 3: Notice that if Condition (4) is replaced by 
E|a;|^ | < L, where i is a fixed positive constant, then we have 

P(Ama.(S„) > L{b + e)) = o(n-*). (94) 
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We only need to consider the matrix t and replace Xij and the above identity, we have 
by b~^/'^Xij to verify the six conditions. 



B. First step truncation and renormalization 

We consider S„T„, where the assumptions of Theorem 3 
are met, except the T„ are assumed nonrandom. Here, to be 
consistent with Chapter 6 of [16], we replace c„ with y„ and 
c with y. 

Notice, from the identity 



max|A^(S„T„)- A^:(S„T„)| 




(103) 
(104) 

(105) 



(95) 



< 



^E|4|E|4|7(|a;i,-| > i^nM (106) 



valid for any nonnegative random variable Y, that (ii) implies 

the fourth moments of the .Xi, exist. 



< 



0. 



(107) 



We will need the following identity later on. For nonnegative 
Y having finite fourth moment, since EF^7(F > y) > 

V^P{Y > y), we have y^P{Y > y) ^ as y ^ oo. Thus, s„ = iwwj Then, by Theorem A. 46 of [16], " 
using integration by parts, we have for any a > " 



3. Reseating. Define Wij — Zij/uij, W„ = {wij)pxn and 



EY'^I{Y > a, 

= a^P(y > o) + lim (-/P(y >y)+ r 4a;^P(F > x)dx) 
y^°° Ja 

(96) 



max|A|(S„T„)- A|(S„T„)| 

k<p 



< ||T^||||n-^(Z, 
1 

7^ 



< 



/■OO 

^P{Y>a)+ / Ax^P{Y>x)dx. 

J a 



w 



0, a.s. 



(108) 
(109) 



(97) because of (94) and the fact that 



We choose 4 such that ?7„\/n t oo, liminf„ rj^\/n > 
0, and 



max 1 1 — A, 



< max[E\xUl{\xij \ > r]y/n) + {E\xij\I{\xij \ > rjy^)) 



(110) 
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< oo, 



(98) 



fe=i 



< 2i;-^{r]y^) maxE\xf M{\xij\) -S- 



where fjk = ri2k . 

1. Truneation. Define yij = Xijl{\xij\ < T]y/n), Y„ 
{yij)pxn and S„ = ^Y„Yi^. We have 

P(X„ ^ Y„, i.o.) 

oo / 2'= + ! 

U U {\x^j\>VnM 

k=m \n='2^+l i^Ptj^t^ 

r,k + l 



which implies that 



maxEi\wij — Zij\ = max 



'if (1 + Aij)2 



(111) 
(112) 



(113) 



(99) 



k=m 



in=2^-\-l i<'2ynj<n 



u 



(100) 

{\xij\ > fjk2''/^} ] (101) 



C. Second truncation and normalization 

We may assume now that the Xij satisfy the six conditions 
of Section B-A with Condition (4) strengthened to E|a;? | = 1 
for all 

Define yij = Xijl{\xij\ < C) — Exijl{\xij\ < C) for some 
large constant C and define Y„ = (?y,:,)pxn, S„ = ^Yn^n- 
Then by Theorem A. 46 of [16], we have 



= lim V P 

m^oo ^ — ^ 
k=m 

c 

= d>yK lim y 22'^P (\X\ > fik2^''^\ = 0. (102) 



max|A^(S„T„)-A;^ (S„T„ 



< IIT^ 



n 



-V2(x„ - Y„) 



< 



k=m 



n 



-1/2CX -Y 



(114) 
(115) 



Since E\x^j - yij\^ < E|4|J(|xy| > C) < M/^(C), this 
2. Centralization. Define Zij = yij — Ej/^ and Z„ = can be made arbitrarily small by making C sufficiently large. 
{zij)pxn and S„ = ^Z„Z5^. Then by Theorem A. 46 of [16] We can then apply (94). 
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The rescaling is the same as given in last Section. Now we 
have 

max|l - afA < 2V'"^(C) maxE|x? |V(^i7), (116) 

ij ij 

which can be made arbitrarily small by making C sufficiently 
large. 

D. Extension of [27] and Chapter 6 of [16] 

The result in [27] on the smallest eigenvalue, Aniin(S„), 
of S„ can be extended with only Assumptions 1), 2), 3) of 
Theorem 3. Indeed, using the two step truncations, we may 
assume the Xij are bounded, with mean and variance 1. 
Following the same steps as in [27], one may prove that when 
y < 1 

Amin(S„)^(l-Vy)', (117) 

almost surely. 

We proceed now to the necessary changes in Chapter 6 in 
[16]. We may now assume the same conditions as in Section 
6.2.1 of [16] on the Xij (except they need not be identically 
distributed), and the bounds appearing there. The changes are 
needed wherever identical distribution was exploited. 

We begin with Page 139, below (6.2.34). We change the 
definition of 6„ to 

1 



" l + n-iEtr(T„D-i)' 
and introduce the quantities 

" l + n-iEtr(T„D-i)' 



(118) 



(119) 



The argument below (6.2.35) is specific for j = 1, but easily 
extends for any j. Following the argument below (6.2.36), we 
can no longer assume E/3i = —zEjS^, nor is bounded, but we 
have 



sup 

uG[a,6] 



n. 



< K. 



We further have 

i>nk = 0k+ Pkbnklk- 

Then, using (6.2.36) 



(120) 



(121) 



J2{bnk - E/3fe) 



<Kn-'Y.v-\E\-fk\^)i (122) 



(123) 



Since bnj - 6„ = 6„6„jE(i trT(D-i - DT^)), we have, 
using Lemma 6.9 of [16] 



n ^ — ^ 



^{bn - b„k) 



<\z\' — , (124) 



and 



bnk - bn\ < K- 



Thus we have 



max sup \bnj\ < K. 
i ue[o,6] 



(125) 
(126) 



For the rest of Section 6.2.3, 6„ is mentioned twice. We need 
to replace it with bnj and the arguments go through without 
any further changes. 

For Section 6.2.4, (6.2.42) needs to be replaced by 



y- 



+ zynE{s„{z)) 



+ tEs„ 
1 " 

= -J2m [r^D^i(Es„T„ + 



(127) 



— Etr(Es„T„+I)-'T„D-i 
n 



(128) 



For the rest of the section, replace subscript 1 with subscript 
k, subscript 2 with subscript j, replace 6i„ with 



bki 



1 



l + n-iEtr(T„Dfej) 



(129) 



and all appearances of subscripts kj assume k ^ j. Fnkj has 
the obvious definition. Replace the summations for j ranging 
from 2 to n with j ^ k. All the bounds derived for k ~ \ 
are true for all k. So we conclude the left side of (6.2.42) is 
bounded by kn~^. 

The rest of Chapter 6 follows without any changes. 
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